##########################################################################################

library(data.table)
library(optparse)
library(Seurat)
library(ArchR)
library(ggthemes)
library(dplyr)

##########################################################################################
option_list <- list(
    make_option(c("--comine_data_all_file"), type = "character"),
    make_option(c("--report_tf_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 所有的细胞的,计算maxdelt
    comine_data_all_file <- "~/20231121_singleMuti/results/qc_atac/testis_combined_peak.combineRNA.qc.Rdata"

    ## 认为关键的TF
    report_tf_file <- "~/20231121_singleMuti/config/TF_upregulated_during_spermatogenesis.csv"

    ## 输出
    out_path <- paste0("~/20231121_singleMuti/results/report_tf")

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

comine_data_all_file <- opt$comine_data_all_file
report_tf_file <- opt$report_tf_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################
## 导入数据

b <- load(comine_data_all_file)
## testis_combined_peak_combineRNA

report_tf <- unique(data.frame(fread(report_tf_file)))
report_tf <- report_tf %>%
group_by( gene ) %>%
summarise( cell_type = paste0(cell_type , collapse = "|") , cluster = paste0(cluster , collapse = "|"))

###########################################################################################
## 转录因子 （TF） 足迹允许预测 TF 在特定位点的精确结合位置。
## 这是因为直接与 TF 结合的 DNA 碱基实际上受到保护，不会转座，而紧邻 TF 结合的 DNA 碱基是可访问的。

projHeme5 <- testis_combined_peak_combineRNA
projHeme5 <- addMotifAnnotations(ArchRProj = projHeme5, motifSet = "cisbp", name = "Motif")
projHeme5 <- addBgdPeaks(projHeme5)
projHeme5 <- addDeviationsMatrix(
  ArchRProj = projHeme5, 
  peakAnnotation = "Motif",
  force = TRUE
)

motifPositions <- getPositions(projHeme5)